home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / MATH1 / ROMB3.PAS < prev    next >
Pascal/Delphi Source File  |  1985-04-03  |  2KB  |  96 lines
  1. program romb3;        { -> 287 }
  2. { integration by the romberg method }
  3.  
  4. const    tol        = 1.0E-4;
  5. var    done        : boolean;
  6.     sumt        : real;
  7.     sum,upper,lower    : real;
  8.  
  9. external procedure cls;
  10.  
  11. function fx(x: real): real;
  12. { find f(x)= 1/sqrt(x); watch out for x=0 }
  13. begin
  14.   fx:=1.0/sqrt(x)
  15. end;
  16.  
  17. procedure romb(function f(x: real): real;
  18.         lower,upper,tol: real;
  19.             var ans: real);
  20. { numerical integration by the Romberg method }
  21. var
  22.     nx            : array[1..16] of integer;
  23.     t            : array[1..136] of real;
  24.     done,error        : boolean;
  25.     pieces,nt,i,ii,n,nn,
  26.     l,ntra,k,m,j        : integer ;
  27.     delta_x,c,sum,fotom,x    : real ;
  28. begin
  29.   done:=false;
  30.   error:=false;
  31.   pieces:=1;
  32.   nx[1]:=1;
  33.   delta_x:=(upper-lower)/pieces;
  34.   c:=(f(lower)+f(upper))*0.5;
  35.   t[1]:=delta_x*c;
  36.   n:=1;
  37.   nn:=2;
  38.   sum:=c;
  39.   repeat
  40.     n:=n+1;
  41.     fotom:=4.0;
  42.     nx[n]:=nn;
  43.     pieces:=pieces*2;
  44.     l:=pieces-1;
  45.     delta_x:=(upper-lower)/pieces;
  46.     { compute trapezoidal sum for 2^(n-1)+1 points }
  47.     for ii:=1 to (l+1) div 2 do
  48.       begin
  49.     i:=ii*2-1;
  50.     x:=lower+i*delta_x;
  51.     sum:=sum+f(x)
  52.       end;
  53.     t[nn]:=delta_x*sum;
  54.     ntra:=nx[n-1];
  55.     k:=n-1;
  56.     { compute n-th row of T array }
  57.     for m:=1 to k do
  58.       begin
  59.     j:=nn+m;
  60.     nt:=nx[n-1]+m-1;
  61.     t[j]:=(fotom*t[j-1]-t[nt])/(fotom-1.0);
  62.     fotom:=fotom*4.0
  63.       end;
  64.     if n>4 then
  65.       begin
  66.     if t[nn+1]<>0.0 then
  67.       if (abs(t[ntra+1]-t[nn+1])<=abs(t[nn+1]*tol))
  68.         or (abs(t[nn-1]-t[j])<=abs(t[j]*tol)) then
  69.           done:=true
  70.       else if n>15 then
  71.         begin
  72.           done:=true;
  73.           error:=true
  74.         end
  75.     end;    { if n>4 }
  76.     nn:=j+1
  77.   until done;
  78.   ans:=t[j]
  79. end;        { ROMBERG }
  80.  
  81. begin        { main program }
  82.   cls;
  83.   lower:=0.1;
  84.   upper:=1.0;
  85.   writeln;
  86.   sumt:=0.0;
  87.   writeln('new area    total area    lower    upper limits');
  88.   repeat
  89.     romb(fx,lower,upper,tol,sum);
  90.     upper:=lower;
  91.     lower:=0.1*upper;
  92.     sumt:=sumt+sum;
  93.     writeln(sum:9:6,'    ',sumt:9:5,'    ',lower,'    ',upper)
  94.   until abs(sum)<tol
  95. end.
  96.